library(tidyverse)
library(skimr)
library(rethinking)

Notes

Introduction

Linear regression describes a measurements mean and variance as the addition of other measurements. It assumes the errors in the primary measurement are of the gaussian distribution.

Normality

1000 individuals conduct a random walk on either side of the football field starting at the 50 yard line (value = 0). This signifies the normal distribution.

INDIVIDUALS <- 1000
TRIALS <- 100

purrr::map(1:INDIVIDUALS, ~ runif(TRIALS, -1, 1)) %>% 
  purrr::map(cumsum) %>% 
  reduce(rbind) %>% 
  as_tibble %>% 
  `colnames<-`(paste0(1:TRIALS, "t")) %>% 
  mutate(individual = row_number()) %>% 
  reshape2::melt("individual") %>% 
  arrange(individual, variable) %>% 
  ggplot(aes(x=variable, y=value, group=individual)) + 
    geom_line(alpha = I(1/20), size=1) +
    geom_boxplot(aes(x=variable, y = value), inherit.aes = FALSE, colour = "yellow", alpha= I(1/40))

The above example shows us that normality is additive.

# Generate 12 numbers between 1 and 1.1 then add them in to one sum
# Do this 1000 times and collect the numbers in a vector
growth <- replicate(1000, sum(runif(12,0,.1)))
# Plot the distribution of thes enumbers. Note that it is normal.
qplot(growth)

We can also show that normality is multiplicative.

# Generate 12 numbers between 1 and 1.1 then multiply them in to one product
# Do this 1000 times and collect the numbers in a vector
growth <- replicate(1000, prod(1 + runif(12,0,.1)))
#plot the distribution of these numbers. Note that it is normal.
qplot(INDIVIDUALS)

Things start to skew at very large deviates, though this can be corrected by log scaling.

big <- replicate(10000, prod(1 + runif(12, 0, .5)))
small <- replicate(10000, prod(1 + runif(12, 0, 0.01)))

data_frame(big = big, 
           small = small, 
           log_big = log(big)) %>% 
  gather %>% 
  ggplot(aes(value)) + 
    geom_histogram() + 
    facet_wrap("key", scales = "free")

A language for describing models

Elements:

An outcome variable: The set of measurements we hope to predict or understand
Likeliehood distribution: For each of these outcome variables, we define the plausibility of that observations.
Predictor variables: Set of other measurements that we hope to use to predict or understand the outcome.

Approach:

Relate the exact shape of the likelihood distribution to the predictor variables via parameters
Choose priors for all the parameters in the model.

Example: A gaussian model of height

We’ll build a regression to predict a human’s hight, which we will model as a Gaussian distribution with two parameters, mean and standard deviation. There are an infinite number of possible Gaussian distributions. We want our Bayesian machine to consider every possible distribution, each defined by a combination of meand and sd, and rank them by posterior plausibility given the data.

data(Howell1)
d <- Howell1
d2 <- d %>% 
  filter(age >= 18)
d2 %>% 
  skim

Skim summary statistics
 n obs: 352 
 n variables: 4 

── Variable type:integer ────────────────────────────────────────────────────────────
 variable missing complete   n mean  sd p0 p25 p50 p75 p100     hist
     male       0      352 352 0.47 0.5  0   0   0   1    1 ▇▁▁▁▁▁▁▇

── Variable type:numeric ────────────────────────────────────────────────────────────
 variable missing complete   n   mean    sd     p0    p25    p50    p75   p100
      age       0      352 352  41.14 15.97  18     28     39     51     88   
   height       0      352 352 154.6   7.74 136.53 148.59 154.31 160.66 179.07
   weight       0      352 352  44.99  6.46  31.07  40.26  44.79  49.29  62.99
     hist
 ▇▇▇▅▂▃▁▁
 ▁▅▇▆▆▃▁▁
 ▂▆▇▇▇▃▂▁

Lets specifiy our model like so. The first line is the likelihood of our outcome, height. The second two lines are priors for the parameters in our likelihood. \[ h_i \sim Normal(\mu, \sigma) \\ \mu \sim Normal(178, 20) \\ \sigma \sim Uniform(0, 50) \]

As height is a physical phenomenon, we can use our intuition to assign a prior.

data_frame(x=100:250, 
           y=purrr::map_dbl(100:250, ~ dnorm(.x, 178, 20))) %>% 
  ggplot(aes(x,y)) + geom_path()

We give the standard deviation a truly flat prior. We really just want to make sure its positive.

data_frame(x=-10:60, 
           y=purrr::map_dbl(-10:60, ~ dunif(.x, 0, 50))) %>% 
  ggplot(aes(x,y)) + geom_path()

Specifying a prior for the parameters of our outcome implicitly give us a prior for our outcome.

sample_mu <- rnorm(1e4, 178, 20)
sample_sigma <- runif(1e4, 0, 50)
prior_h <- rnorm(1e4, sample_mu, sample_sigma)
qplot(prior_h)

Lets us grid approximation to generate the posterior. To do so, we enumerate a list of candidate values for our parameters and store them in post. We then calculate the log likelihood for each of those parameters given all the heights we have seen and store them in the column LL. We use log transformation to avoid rounding errors. Then we multiply by the priors (log addition is the same as multiplying)

# Generate candidate parameters
mu_list <- seq(from=140, to=160, length.out=200)
sigma_list <- seq(from=4, to=9, length.out=200)
post <- expand.grid(mu=mu_list, sigma=sigma_list)

# Calculate the likelihood of every possible
# combination candidate parameters
post$LL <- sapply(1:nrow(post), function(i) {
  sum(dnorm(
    d2$height,
    mean=post$mu[i],
    sd=post$sigma[i],
    log=TRUE))
  })

# Multiply priors
post$prod <- post$LL + dnorm(post$mu, 178, 20, TRUE) +
  dunif(post$sigma, 0, 50, TRUE)

# Calculated scaled posterior
post$prob <- exp(post$prod - max(post$prod))

# Contour plot of probabilities
ggplot(post, aes(x=mu, y=sigma, z=prob)) + 
  geom_contour(aes(colour=stat(level))) + 
  coord_cartesian(xlim=c(140, 160), ylim=c(4, 9))

Now lets sample from the posterior.

# Sample candidate values by posterior
sample_rows <- sample(1:nrow(post), size=1e4, replace=TRUE, prob=post$prob)

# Plot their 2d histogram
data_frame(
  sample_mu = post$mu[sample_rows],
  sample_sigma = post$sigma[sample_rows]
) %>% 
  ggplot(aes(sample_mu, sample_sigma)) +
    geom_bin2d(bins=30) + 
    coord_cartesian(xlim=c(140, 160), ylim=c(4, 9))

Lets also examine the margins.

require(gridExtra)

Loading required package: gridExtra

Attaching package: ‘gridExtra’

The following object is masked from ‘package:dplyr’:

    combine

p1 <- post$mu[sample_rows] %>% 
  qplot + ggtitle("mu")
p2 <- post$sigma[sample_rows] %>% 
  qplot + ggtitle("sigma")

grid.arrange(p1, p2, ncol=2)

map from the rethinking package allows us to solve for the posterior using quadratic approximation instead of grid approximations. It shows us gaussian approximations for each parameter’s marginal distribution.

flist <- alist(
  height ~ dnorm(mu, sigma),
  mu ~ dnorm(178, 20),
  sigma ~ dunif(0, 50)
)
m4.1 <- rethinking::map(flist, data=d2)
precis(m4.1)

Variance covariance matrix tells us how each parameter relates to every other parameter in the posterior distribution.

vcov(m4.1)

                mu        sigma
mu    0.1697392273 0.0002182711
sigma 0.0002182711 0.0849053439

# This can be factored in to a vector of variances
diag(vcov(m4.1))

        mu      sigma 
0.16973923 0.08490534

# Or a correlation matrix
cov2cor(vcov(m4.1))

               mu       sigma
mu    1.000000000 0.001818183
sigma 0.001818183 1.000000000

Lets draw from the posterior by sampling vectors of values from a multi-dimensional gaussian distribution.

post <- extract.samples(m4.1, n=1e4)
precis(post)

quap posterior: 10000 samples from m4.1

Adding a predictor variable, weight

It looks like weight will be a good predictor for height in adults.

ggplot(d2, aes(x=weight, y=height)) + 
  geom_point(shape=1) +
  geom_smooth(method=lm)

The goal here is to make the parameter for the mean of a Gaussian distribution, \(\mu\), into a linear function of the predictor variable. This encodes the assumption that the predictor variable has a perfectly constant and additive relationship to the mean of the outcome.

m4.3 <- rethinking::map(
  alist(
    height ~ dnorm(mu, sigma),
    mu <- a + b*weight,
    a ~ dnorm(156, 100),
    b ~ dnorm(0, 10),
    sigma ~ dunif(0, 50)),
  data = d2)
# Display table of estimates
precis(m4.3, corr=TRUE)

The table above tells us that a person 1 kg heavier is expected to be .9 cm taller. It also tells us that a person with 0 weight is 114 cm tall, which is of course nonsense; This is why it is usually important to have very weak priors for intercepts in many cases. Finally, it explains that 95% of plausible heights lie within 10 cm of the mean.

Centering is the procedure of subtracting the mean of a variable from each value.

d2$weight.c <- d2$weight - mean(d2$weight)

m4.4 <- rethinking::map(
  alist(
    height ~ dnorm(mu, sigma),
    mu <- a + b*weight.c,
    a ~ dnorm(178, 100),
    b ~ dnorm(0, 10),
    sigma ~ dunif(0, 50)),
  data = d2)
# Display table of estimates
precis(m4.4, corr=TRUE)

By centering our predictor variable, the interpretation of the intercept has now become the expected value of the outcome when the predictor is at its average value. This makes interpreting the intercept a lot easier.

Visualizing our estimates

We can plot our maximum a posteriori fit here

p <- ggplot(d2, aes(x=weight, y=height)) +
  geom_point(shape=1, size=2) +
  geom_abline(aes(intercept = coef(m4.3)["a"],
                  slope = coef(m4.3)["b"]))
p

This line only represents one piece of the posterior. We could instead sample values of \(\alpha\) and \(\beta\) from the posterior to generate many lines to build bounds of uncertainty. We’ll do this with a much smaller dataset to magnify effects.

N <- 10
dN <- d2[ 1:N , ]
mN <- rethinking::map(
  alist(
    height ~ dnorm( mu , sigma ) ,
    mu <- a + b*weight ,
    a ~ dnorm( 178 , 100 ) ,
    b ~ dnorm( 0 , 10 ) ,
    sigma ~ dunif( 0 , 50 )
) , data=dN )

post <- extract.samples(mN, n=20)
head(post)

p <- ggplot(dN, aes(weight, height)) +
  geom_point(shape=1)

for (i in 1:20) {
  p <- p + geom_abline(intercept = post$a[i], slope = post$b[i], alpha = I(3/10))
}

p

To do this more generally:

Sample parameter values from the posterior.
Use the parameter values to create predictions of the parameters of interest, including the final outcome.
Use summary functions like mean, HDOP, PI to find averages and lower and upper bounds of our uncertainty.

# Draw parameter values from the posterior distribution
posterior_samples <- extract.samples(m4.3)

# Create out x-axis of interest
weights_axis <- 25:70

# Function to calculate fit values for \mu
mu_link <- function(weight) post$a + post$b*weight
mu <- purrr::map(weights_axis, mu_link)

# Function to calculate fit heights
# This incorporates variability as well
simulate <- function(weight, parameter_samples) {
  rnorm(n = nrow(parameter_samples),
        mean = parameter_samples$a + parameter_samples$b * weight,
        sd = parameter_samples$sigma)}
simulated_heights <- purrr::map(weights_axis, ~ simulate(.x, posterior_samples))

# Combine all this information to plot
data_frame(
  # X axis
  weight = weights_axis,
  # Estimates for mean height
  mu_mean = map_dbl(mu, mean),
  mu_lower = map_dbl(mu, ~HPDI(.x, .89)[1]),
  mu_upper = map_dbl(mu, ~HPDI(.x, .89)[2]),
  # Estimates for height
  height_mean = map_dbl(simulated_heights, mean),
  height_low = map_dbl(simulated_heights, ~HPDI(.x, .89)[1]),
  height_high = map_dbl(simulated_heights, ~HPDI(.x, .89)[2])) %>% 
  ggplot(data=.) +
  geom_point(inherit.aes = FALSE, data = d2,
             aes(x=weight, y=height), shape=1) +
  geom_line(aes(x=weight, y=mu_mean)) +
  geom_ribbon(aes(x = weight, ymin=mu_lower, ymax=mu_upper), alpha=.5) +
  geom_ribbon(aes(x = weight, ymin = height_low, ymax = height_high), alpha=.3) +
  ggtitle("The effect of weight on height", subtitle = "Using 89% highest posterior density intervals") +
  coord_cartesian(xlim=c(30, 60),ylim=c(135, 180))

Polynomial Regression

Standardizing our predictors means scaling but also dividing by the standard deviation. This changes our interpretation of continuous variables to mean the unit change in outcome from a cahnge of one standard deviation.

# Standardize continuous variable
standardize <- function(vec) {
  centered <- vec - mean(vec)
  centered / sd(vec)
}
d$weight.s <- standardize(d$weight)

# Define another polynomial term
d$weight.s2 <- d$weight.s^2

m4.5 <- map(
  alist(
    # Outcome
    height ~ dnorm(mu, sigma),
    # Regression specification
    mu <- a + b1*weight.s + b2*weight.s2,
    a ~ dnorm(178 , 100),
    # These priors are very weak!
    b1 ~ dnorm(0 , 10),
    b2 ~ dnorm(0 , 10),
    sigma ~ dunif(0 , 50)),
  data=d )
precis( m4.5 )

B-Splines

Splines are an alternative to modeling non-linearity. We’ll show this with some cherry blossom data.

library(rethinking)
library(rcartocolor)
colour_theme <- "BurgYl"
palette <- carto_pal(7, colour_theme)

data(cherry_blossoms)
d <- cherry_blossoms

d %>% 
  select(year, temp) %>% 
  ggplot(aes(year, temp)) +
    geom_line(colour = palette[7]) +
    theme_burgyl()

NA

Splines work by defining knot points and fitting local models between them. A common way to pick these knot points is to evenly space them across the probability density of the variable. Cross validation can ultimately be used to validate the number and spacing.

d2 <- d %>% 
  filter(!is.na(temp))
num_knots <- 15
knot_list <- quantile(d2$year, probs = seq(0, 1, length.out = num_knots))
d2 %>% 
  ggplot(aes(year)) +
  geom_density(colour = "transparent", fill = palette[5]) +
  geom_point(data = data_frame(y= rep(0, length(knot_list)), x=knot_list),
             mapping = aes(x, y), 
             size = 3, colour = palette[7]) +
  theme_burgyl() + ggtitle("Knot point allocation")

Now we pick the polynomial degree. This determines how many basis parameters at once are interacting with a given point. The splines library in R can help us with this. The code below will create 17 columns in place of year that result from its change of basis.

posterior_summary(fit.spline)

                 Estimate   Est.Error         Q2.5        Q97.5
b_Intercept    6.32232353 0.239223616    5.8391272    6.7829623
b_year_1       0.09842464 0.260051555   -0.4001889    0.6075442
b_year_2       0.24137556 0.275882205   -0.3072947    0.7821862
b_year_3       1.19993220 0.269603691    0.6743946    1.7293925
b_year_4      -0.82019659 0.255732821   -1.3066008   -0.3167264
b_year_5       0.10171999 0.254785921   -0.4080378    0.6020521
b_year_6      -1.40488733 0.252521101   -1.8872237   -0.8965606
b_year_7       1.15627188 0.255715220    0.6617547    1.6692874
b_year_8      -1.91663644 0.250916348   -2.4033165   -1.4195968
b_year_9       2.34684446 0.254658438    1.8454632    2.8576482
b_year_10     -2.31889187 0.252143992   -2.8041959   -1.8086396
b_year_11      0.93339356 0.254833356    0.4261467    1.4416483
b_year_12     -1.61366922 0.254398401   -2.1167049   -1.1142676
b_year_13      0.17894803 0.253669626   -0.3172705    0.6873232
b_year_14     -1.24132292 0.256076052   -1.7392509   -0.7453285
b_year_15      0.03278374 0.268382189   -0.4885031    0.5606930
b_year_16      0.99338364 0.270930809    0.4639197    1.5279240
b_year_17      2.03762622 0.265267132    1.5035638    2.5580490
sigma          0.34736003 0.007543462    0.3330691    0.3622407
lp__        -441.28700233 3.153867138 -448.2382129 -436.2200223

Extract the predicted means and the error term.

mu = as_data_frame(fitted(fit.spline))
temp_hat = as_data_frame(predict(fit.spline))

(dimensions <- list(mu = dim(mu), 
                   temp_hat = dim(temp_hat)))

$mu
[1] 1124    4

$temp_hat
[1] 1124    4

Plot the posterior predictions.

data_frame(
  temp = filter(d, !is.na(temp))$temp,
  year = filter(d, !is.na(temp))$year,
  temp_mu = mu$Estimate,
  temp_mu_lb = mu$Q2.5,
  temp_mu_ub = mu$Q97.5,
  temp_hat_p = temp_hat$Estimate,
  temp_hat_lb = temp_hat$Q2.5,
  temp_hat_ub = temp_hat$Q97.5) %>% 
  ggplot(aes(x = year, y = temp)) +
    geom_jitter(alpha = .4, colour = palette[6]) +
    geom_ribbon(aes(ymin = temp_mu_lb, ymax = temp_mu_ub), fill = alpha(palette[4], .7)) +
    geom_ribbon(aes(ymin = temp_hat_lb, ymax = temp_hat_ub), fill = alpha(palette[4], .4)) +
    geom_line(aes(y = temp_hat_p), colour = palette[7]) +
    theme_burgyl() +
    ggtitle("Posterior prediction check",
            "Fitting temperature to a 15 knotted cubic spline of year")

Actually, brms can using the s or t2 functions from the mgcv package to accomplish the same thing.

mu = as_data_frame(fitted(fit.spline.2))
temp_hat = as_data_frame(predict(fit.spline.2))

data_frame(
  temp = filter(d, !is.na(temp))$temp,
  year = filter(d, !is.na(temp))$year,
  temp_mu = mu$Estimate,
  temp_mu_lb = mu$Q2.5,
  temp_mu_ub = mu$Q97.5,
  temp_hat_p = temp_hat$Estimate,
  temp_hat_lb = temp_hat$Q2.5,
  temp_hat_ub = temp_hat$Q97.5) %>% 
  ggplot(aes(x = year, y = temp)) +
    geom_jitter(alpha = .4, colour = palette[6]) +
    geom_ribbon(aes(ymin = temp_mu_lb, ymax = temp_mu_ub), fill = alpha(palette[4], .7)) +
    geom_ribbon(aes(ymin = temp_hat_lb, ymax = temp_hat_ub), fill = alpha(palette[4], .4)) +
    geom_line(aes(y = temp_hat_p), colour = palette[7]) +
    theme_burgyl() +
    ggtitle("Posterior prediction check",
            "Fitting temperature to a 15 knotted cubic spline of year")

On transformations

Centering and scaling a continuous variable increases interpretability. Demeaning a predictor variable subtracts every value by it’s mean. By demeaning a predictor variables, the intercept represents the value of the outcome when all variables are at their mean. One standardizes a variable by dividing every value by the variables standard deviation. Standardizing causes the interpretation of coefficients to be in terms of unit changes in SD of the variable.

Log transforming the outcome variable makes y our coefficients represent percent increases in the outcome. Log transforming a predictor variable describes how a single percentage point increase in the predictor affects y. Log transforming both the predictor and outcome describes how percentage changes in one create percentage changes in the other.

Homework

Medium

4M1

For the model definition below, simulate observed heights from the prior.

num_observations <- 1e4

prior_mean <- rnorm(num_observations, 0, 10)
prior_sd <- runif(num_observations, 0, 10)

prior_heights <- rnorm(num_observations, prior_mean, prior_sd)

qplot(prior_heights)

4M2

model <- alist(
  y ~ dnorm(mu, sigma),
  mu ~ dnorm(0, 10),
  sigma ~ dunif(0, 10)
)

Hard

4H1

data("Howell1")
d <- Howell1
d$weight <- standardize(d$weight)
d

model <- map(
  alist(
    height ~ dnorm(mu, sigma),
    mu <- a + b * weight,
    a ~ dnorm(mean = 0, sd = 100),
    b ~ dnorm(mean = 0, sd = 10),
    sigma ~ dunif(min = 0, max = 64)
  ),
  data = d
)
model

# Standardize our input
new_weights <- standardize(c(46.95, 43.72, 64.78, 32.59, 54.63))

# Draw parameter values from the posterior distribution
posterior_samples <- extract.samples(model)

# Function to calculate fit heights
simulate_heights <- function(weight, parameter_samples) {
  rnorm(
    n = nrow(parameter_samples),
    mean = parameter_samples$a + parameter_samples$b * weight,
    sd = parameter_samples$sigma)}

# Generate a prediction for each new weight
# Use every parameter set in the posterior to do so
simulated_heights <- purrr::map(new_weights, ~ simulate(.x, posterior_samples))

data_frame(
  individual = 1:5,
  weight = c(46.95, 43.72, 64.78, 32.59, 54.63),
  expected_height = map_dbl(simulated_heights, mean),
  lower_pi = map_dbl(simulated_heights, ~PI(.x)[1]),
  upper_pi = map_dbl(simulated_heights, ~PI(.x)[2])
)

4H2

data("Howell1")
d <- Howell1
d %>% 
  filter(age < 18) %>% 
  ggplot(aes(x = weight, y = height)) +
  geom_point() + geom_smooth(method=lm) +
  geom_smooth(colour = "orange")

d %>% 
  filter(age < 18) %>% 
  map(alist(
    height ~ dnorm(mu, sigma),
    mu <- a + b * weight,
    a ~ dnorm(mean = 100, sd = 100),
    b ~ dnorm(mean = 0, sd = 10),
    sigma ~ dunif(min = 0, max = 50)
  ),
  data = .) ->
model

precis(model)

A one unit increase in weight increases height by 2.72, therefore a 10 unit increase will move height by 27.2.

# Sample the poseterior
posterior_samples <- extract.samples(model)

# Create x-axis
new_weights <- seq(from = 4, to = 45, length.out = 100)

# Calculate the mean
fn_mean <- function(weight) {posterior_samples$a + posterior_samples$b * weight}
fit_means <- purrr::map(new_weights, fn_mean)

# Calculate heights
fn_height <- function(weight) {
  rnorm(
    n = nrow(posterior_samples),
    mean = posterior_samples$a + posterior_samples$b * weight,
    sd = posterior_samples$sigma)
}
fit_heights <- purrr::map(new_weights, fn_height)

data_frame(
  weights = new_weights,
  means = map_dbl(fit_means, mean),
  means_low = map_dbl(fit_means, ~ HPDI(.x, .89)[1]),
  means_high = map_dbl(fit_means, ~ HPDI(.x, .89)[2]),
  heights = map_dbl(fit_heights, mean),
  heights_low = map_dbl(fit_heights, ~ HPDI(.x, .89)[1]),
  heights_high = map_dbl(fit_heights, ~ HPDI(.x, .89)[2]),
) %>% 
ggplot(aes(x = weights)) +
  geom_line(aes(y = means)) +
  geom_ribbon(aes(ymin = means_low, ymax = means_high), alpha = .4) +
  geom_ribbon(aes(ymin = heights_low, ymax = heights_high), alpha = .2) +
  geom_point(inherit.aes = FALSE, data = filter(d, age < 18),
             aes(x=weight, y = height), shape = 1) +
  ylab("height") + xlab("weight") + theme_bw() +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())

BRMS

data(Howell1)
d <- Howell1
d2 <- d %>% 
  filter(age >= 18)
d2 %>% 
  skim

library(brms)

b4.1 <-
  brm(data = d2, family = gaussian,
      height ~ 1,
      prior = c(prior(normal(178, 20), class = Intercept),
                prior(cauchy(0, 1), class = sigma)),
      iter = 31000, warmup = 30000, chains = 4, cores = 4)
plot(b4.1)

Extract variance covariance matrix.

# Extract samples from the posterior
post <- posterior_samples(b4.1)
# Covariance matrix
cov(post[, 1:2])

Lets extract the HDPI

library(tidybayes)

post %>% 
  select(-lp__) %>% 
  gather(parameter) %>% 
  group_by(parameter) %>% 
  median_hdi(value)

Now lets add a predictor.

d3 <- d2 %>% 
  mutate(weight = (weight - mean(weight))/sd(weight)) %>% 
  select(height, weight)

b4.3 <- 
  brm(data = d3, family = gaussian,
      height ~ 1 + weight,
      prior = c(prior(normal(156, 100), class = Intercept),
                prior(normal(0, 10), class = b),
                prior(uniform(0, 50), class = sigma)),
      iter = 46000, warmup = 45000, chains = 4, cores = 4,
      control = list(adapt_delta = 0.8, 
                     max_treedepth = 10))

plot(b4.3)

Examine covariance

pairs(b4.3)

Creating predictions

post <- posterior_samples(b4.3)

mu_at_50 <-
  post %>% 
  transmute(mu_at_50 = b_Intercept + b_weight * 50)

mu_at_50 %>%
  ggplot(aes(x = mu_at_50)) +
  geom_density(size = 0, fill = "grey75") +
  stat_pointintervalh(aes(y = 0), 
                      point_interval = mode_hdi, .width = .95) +
  scale_y_continuous(NULL, breaks = NULL) +
  labs(x = expression(mu["height | weight = 50"])) +
  theme_classic()

Now we extract complete predictions including the mean and added variance.

weight_seq <- tibble(weight = seq(from = 25, to = 70, by = 1))

pred_height <-
  predict(b4.3,
          newdata = weight_seq) %>%
  as_tibble() %>%
  bind_cols(weight_seq)
  
pred_height %>%
  slice(1:6)

ggplot(pred_height, aes(weight, Estimate)) +
  geom_ribbon(aes(ymin = Q2.5, ymax = Q97.5), fill = "grey75") +
  geom_line() +
  theme(panel.grid = element_blank())

We can overlay this with the mean interval too

mu_summary <-
  fitted(b4.3, 
         newdata = weight_seq) %>%
  as_tibble() %>%
  bind_cols(weight_seq)

ggplot(mu_summary, aes(weight, Estimate)) +
  geom_ribbon(aes(ymin = Q2.5, ymax = Q97.5), fill = "grey75") +
  geom_line() +
  theme(panel.grid = element_blank())

inner_join(pred_height, mu_summary, by = "weight", suffix=c("_height", "_mean")) %>% str

inner_join(pred_height, mu_summary, by = "weight", suffix=c("_height", "_mean")) %>% 
  ggplot(aes(x = weight, y = Estimate_height)) +
  geom_ribbon(aes(ymin = Q2.5_height, ymax = Q97.5_height), fill = "grey25") +
  geom_ribbon(aes(ymin = Q2.5_mean, ymax = Q97.5_mean), fill = "grey45") +
  geom_line() +
  theme(panel.grid = element_blank())

LS0tCnRpdGxlOiAiQ2hhcHRlciA0IC0gTGluZWFyIE1vZGVscyIKb3V0cHV0OiBodG1sX25vdGVib29rCi0tLQoKYGBge3IgbWVzc2FnZT1GQUxTRSwgd2FybmluZz1GQUxTRX0KbGlicmFyeSh0aWR5dmVyc2UpCmxpYnJhcnkoc2tpbXIpCmxpYnJhcnkocmV0aGlua2luZykKYGBgCgojIE5vdGVzCgojIyBJbnRyb2R1Y3Rpb24KTGluZWFyIHJlZ3Jlc3Npb24gZGVzY3JpYmVzIGEgbWVhc3VyZW1lbnRzIG1lYW4gYW5kIHZhcmlhbmNlIGFzIHRoZSBhZGRpdGlvbiBvZiBvdGhlciBtZWFzdXJlbWVudHMuIEl0IGFzc3VtZXMgdGhlIGVycm9ycyBpbiB0aGUgcHJpbWFyeSBtZWFzdXJlbWVudCBhcmUgb2YgdGhlIGdhdXNzaWFuIGRpc3RyaWJ1dGlvbi4KCiMjIE5vcm1hbGl0eQoKMTAwMCBpbmRpdmlkdWFscyBjb25kdWN0IGEgcmFuZG9tIHdhbGsgb24gZWl0aGVyIHNpZGUgb2YgdGhlIGZvb3RiYWxsIGZpZWxkIHN0YXJ0aW5nIGF0IHRoZSA1MCB5YXJkIGxpbmUgKHZhbHVlID0gMCkuIFRoaXMgc2lnbmlmaWVzIHRoZSBub3JtYWwgZGlzdHJpYnV0aW9uLgpgYGB7ciBmaWcud2lkdGg9MTR9CklORElWSURVQUxTIDwtIDEwMDAKVFJJQUxTIDwtIDEwMAoKcHVycnI6Om1hcCgxOklORElWSURVQUxTLCB+IHJ1bmlmKFRSSUFMUywgLTEsIDEpKSAlPiUgCiAgcHVycnI6Om1hcChjdW1zdW0pICU+JSAKICByZWR1Y2UocmJpbmQpICU+JSAKICBhc190aWJibGUgJT4lIAogIGBjb2xuYW1lczwtYChwYXN0ZTAoMTpUUklBTFMsICJ0IikpICU+JSAKICBtdXRhdGUoaW5kaXZpZHVhbCA9IHJvd19udW1iZXIoKSkgJT4lIAogIHJlc2hhcGUyOjptZWx0KCJpbmRpdmlkdWFsIikgJT4lIAogIGFycmFuZ2UoaW5kaXZpZHVhbCwgdmFyaWFibGUpICU+JSAKICBnZ3Bsb3QoYWVzKHg9dmFyaWFibGUsIHk9dmFsdWUsIGdyb3VwPWluZGl2aWR1YWwpKSArIAogICAgZ2VvbV9saW5lKGFscGhhID0gSSgxLzIwKSwgc2l6ZT0xKSArCiAgICBnZW9tX2JveHBsb3QoYWVzKHg9dmFyaWFibGUsIHkgPSB2YWx1ZSksIGluaGVyaXQuYWVzID0gRkFMU0UsIGNvbG91ciA9ICJ5ZWxsb3ciLCBhbHBoYT0gSSgxLzQwKSkKYGBgCgpUaGUgYWJvdmUgZXhhbXBsZSBzaG93cyB1cyB0aGF0IG5vcm1hbGl0eSBpcyBhZGRpdGl2ZS4gCmBgYHtyfQojIEdlbmVyYXRlIDEyIG51bWJlcnMgYmV0d2VlbiAxIGFuZCAxLjEgdGhlbiBhZGQgdGhlbSBpbiB0byBvbmUgc3VtCiMgRG8gdGhpcyAxMDAwIHRpbWVzIGFuZCBjb2xsZWN0IHRoZSBudW1iZXJzIGluIGEgdmVjdG9yCmdyb3d0aCA8LSByZXBsaWNhdGUoMTAwMCwgc3VtKHJ1bmlmKDEyLDAsLjEpKSkKIyBQbG90IHRoZSBkaXN0cmlidXRpb24gb2YgdGhlcyBlbnVtYmVycy4gTm90ZSB0aGF0IGl0IGlzIG5vcm1hbC4KcXBsb3QoZ3Jvd3RoKQpgYGAKCldlIGNhbiBhbHNvIHNob3cgdGhhdCBub3JtYWxpdHkgaXMgbXVsdGlwbGljYXRpdmUuCmBgYHtyfQojIEdlbmVyYXRlIDEyIG51bWJlcnMgYmV0d2VlbiAxIGFuZCAxLjEgdGhlbiBtdWx0aXBseSB0aGVtIGluIHRvIG9uZSBwcm9kdWN0CiMgRG8gdGhpcyAxMDAwIHRpbWVzIGFuZCBjb2xsZWN0IHRoZSBudW1iZXJzIGluIGEgdmVjdG9yCmdyb3d0aCA8LSByZXBsaWNhdGUoMTAwMCwgcHJvZCgxICsgcnVuaWYoMTIsMCwuMSkpKQojcGxvdCB0aGUgZGlzdHJpYnV0aW9uIG9mIHRoZXNlIG51bWJlcnMuIE5vdGUgdGhhdCBpdCBpcyBub3JtYWwuCnFwbG90KElORElWSURVQUxTKQpgYGAKClRoaW5ncyBzdGFydCB0byBza2V3IGF0IHZlcnkgbGFyZ2UgZGV2aWF0ZXMsIHRob3VnaCB0aGlzIGNhbiBiZSBjb3JyZWN0ZWQgYnkgbG9nIHNjYWxpbmcuCmBgYHtyIH0KYmlnIDwtIHJlcGxpY2F0ZSgxMDAwMCwgcHJvZCgxICsgcnVuaWYoMTIsIDAsIC41KSkpCnNtYWxsIDwtIHJlcGxpY2F0ZSgxMDAwMCwgcHJvZCgxICsgcnVuaWYoMTIsIDAsIDAuMDEpKSkKCmRhdGFfZnJhbWUoYmlnID0gYmlnLCAKICAgICAgICAgICBzbWFsbCA9IHNtYWxsLCAKICAgICAgICAgICBsb2dfYmlnID0gbG9nKGJpZykpICU+JSAKICBnYXRoZXIgJT4lIAogIGdncGxvdChhZXModmFsdWUpKSArIAogICAgZ2VvbV9oaXN0b2dyYW0oKSArIAogICAgZmFjZXRfd3JhcCgia2V5Iiwgc2NhbGVzID0gImZyZWUiKQpgYGAKCiMjIEEgbGFuZ3VhZ2UgZm9yIGRlc2NyaWJpbmcgbW9kZWxzCgpFbGVtZW50czoKCjEuICoqQW4gb3V0Y29tZSB2YXJpYWJsZSoqOiBUaGUgc2V0IG9mIG1lYXN1cmVtZW50cyB3ZSBob3BlIHRvIHByZWRpY3Qgb3IgdW5kZXJzdGFuZAoyLiAqKkxpa2VsaWVob29kIGRpc3RyaWJ1dGlvbioqOiBGb3IgZWFjaCBvZiB0aGVzZSBvdXRjb21lIHZhcmlhYmxlcywgd2UgZGVmaW5lIHRoZSBwbGF1c2liaWxpdHkgb2YgdGhhdCBvYnNlcnZhdGlvbnMuCjMuICoqUHJlZGljdG9yIHZhcmlhYmxlcyoqOiBTZXQgb2Ygb3RoZXIgbWVhc3VyZW1lbnRzIHRoYXQgd2UgaG9wZSB0byB1c2UgdG8gcHJlZGljdCBvciB1bmRlcnN0YW5kIHRoZSBvdXRjb21lLgoKQXBwcm9hY2g6CgoxLiBSZWxhdGUgdGhlIGV4YWN0IHNoYXBlIG9mIHRoZSBsaWtlbGlob29kIGRpc3RyaWJ1dGlvbiB0byB0aGUgcHJlZGljdG9yIHZhcmlhYmxlcyB2aWEgcGFyYW1ldGVycwoyLiBDaG9vc2UgcHJpb3JzIGZvciBhbGwgdGhlIHBhcmFtZXRlcnMgaW4gdGhlIG1vZGVsLgoKCiMjIEV4YW1wbGU6IEEgZ2F1c3NpYW4gbW9kZWwgb2YgaGVpZ2h0CgpXZSdsbCBidWlsZCBhIHJlZ3Jlc3Npb24gdG8gcHJlZGljdCBhIGh1bWFuJ3MgaGlnaHQsIHdoaWNoIHdlIHdpbGwgbW9kZWwgYXMgYSBHYXVzc2lhbiBkaXN0cmlidXRpb24gd2l0aCB0d28gcGFyYW1ldGVycywgbWVhbiBhbmQgc3RhbmRhcmQgZGV2aWF0aW9uLiBUaGVyZSBhcmUgYW4gaW5maW5pdGUgbnVtYmVyIG9mIHBvc3NpYmxlIEdhdXNzaWFuIGRpc3RyaWJ1dGlvbnMuIFdlIHdhbnQgb3VyIEJheWVzaWFuIG1hY2hpbmUgdG8gY29uc2lkZXIgZXZlcnkgcG9zc2libGUgZGlzdHJpYnV0aW9uLCBlYWNoIGRlZmluZWQgYnkgYSBjb21iaW5hdGlvbiBvZiBtZWFuZCBhbmQgc2QsIGFuZCByYW5rIHRoZW0gYnkgcG9zdGVyaW9yIHBsYXVzaWJpbGl0eSBnaXZlbiB0aGUgZGF0YS4KCmBgYHtyfQpkYXRhKEhvd2VsbDEpCmQgPC0gSG93ZWxsMQpkMiA8LSBkICU+JSAKICBmaWx0ZXIoYWdlID49IDE4KQpkMiAlPiUgCiAgc2tpbQpgYGAKCkxldHMgc3BlY2lmaXkgb3VyIG1vZGVsIGxpa2Ugc28uIFRoZSBmaXJzdCBsaW5lIGlzIHRoZSBsaWtlbGlob29kIG9mIG91ciBvdXRjb21lLCBoZWlnaHQuIFRoZSBzZWNvbmQgdHdvIGxpbmVzIGFyZSBwcmlvcnMgZm9yIHRoZSBwYXJhbWV0ZXJzIGluIG91ciBsaWtlbGlob29kLgokJApoX2kgXHNpbSBOb3JtYWwoXG11LCBcc2lnbWEpIFxcClxtdSBcc2ltIE5vcm1hbCgxNzgsIDIwKSBcXApcc2lnbWEgXHNpbSBVbmlmb3JtKDAsIDUwKQokJAoKQXMgaGVpZ2h0IGlzIGEgcGh5c2ljYWwgcGhlbm9tZW5vbiwgd2UgY2FuIHVzZSBvdXIgaW50dWl0aW9uIHRvIGFzc2lnbiBhIHByaW9yLgpgYGB7cn0KZGF0YV9mcmFtZSh4PTEwMDoyNTAsIAogICAgICAgICAgIHk9cHVycnI6Om1hcF9kYmwoMTAwOjI1MCwgfiBkbm9ybSgueCwgMTc4LCAyMCkpKSAlPiUgCiAgZ2dwbG90KGFlcyh4LHkpKSArIGdlb21fcGF0aCgpCmBgYAoKV2UgZ2l2ZSB0aGUgc3RhbmRhcmQgZGV2aWF0aW9uIGEgdHJ1bHkgZmxhdCBwcmlvci4gV2UgcmVhbGx5IGp1c3Qgd2FudCB0byBtYWtlIHN1cmUgaXRzIHBvc2l0aXZlLgpgYGB7cn0KZGF0YV9mcmFtZSh4PS0xMDo2MCwgCiAgICAgICAgICAgeT1wdXJycjo6bWFwX2RibCgtMTA6NjAsIH4gZHVuaWYoLngsIDAsIDUwKSkpICU+JSAKICBnZ3Bsb3QoYWVzKHgseSkpICsgZ2VvbV9wYXRoKCkKYGBgCgpTcGVjaWZ5aW5nIGEgcHJpb3IgZm9yIHRoZSBwYXJhbWV0ZXJzIG9mIG91ciBvdXRjb21lIGltcGxpY2l0bHkgZ2l2ZSB1cyBhIHByaW9yIGZvciBvdXIgb3V0Y29tZS4KYGBge3J9CnNhbXBsZV9tdSA8LSBybm9ybSgxZTQsIDE3OCwgMjApCnNhbXBsZV9zaWdtYSA8LSBydW5pZigxZTQsIDAsIDUwKQpwcmlvcl9oIDwtIHJub3JtKDFlNCwgc2FtcGxlX211LCBzYW1wbGVfc2lnbWEpCnFwbG90KHByaW9yX2gpCmBgYAoKTGV0cyB1cyBncmlkIGFwcHJveGltYXRpb24gdG8gZ2VuZXJhdGUgdGhlIHBvc3Rlcmlvci4gVG8gZG8gc28sIHdlIGVudW1lcmF0ZSBhIGxpc3Qgb2YgY2FuZGlkYXRlIHZhbHVlcyBmb3Igb3VyIHBhcmFtZXRlcnMgYW5kIHN0b3JlIHRoZW0gaW4gcG9zdC4gV2UgdGhlbiBjYWxjdWxhdGUgdGhlIGxvZyBsaWtlbGlob29kIGZvciBlYWNoIG9mIHRob3NlIHBhcmFtZXRlcnMgZ2l2ZW4gYWxsIHRoZSBoZWlnaHRzIHdlIGhhdmUgc2VlbiBhbmQgc3RvcmUgdGhlbSBpbiB0aGUgY29sdW1uIExMLiBXZSB1c2UgbG9nIHRyYW5zZm9ybWF0aW9uIHRvIGF2b2lkIHJvdW5kaW5nIGVycm9ycy4gVGhlbiB3ZSBtdWx0aXBseSBieSB0aGUgcHJpb3JzIChsb2cgYWRkaXRpb24gaXMgdGhlIHNhbWUgYXMgbXVsdGlwbHlpbmcpCmBgYHtyfQojIEdlbmVyYXRlIGNhbmRpZGF0ZSBwYXJhbWV0ZXJzCm11X2xpc3QgPC0gc2VxKGZyb209MTQwLCB0bz0xNjAsIGxlbmd0aC5vdXQ9MjAwKQpzaWdtYV9saXN0IDwtIHNlcShmcm9tPTQsIHRvPTksIGxlbmd0aC5vdXQ9MjAwKQpwb3N0IDwtIGV4cGFuZC5ncmlkKG11PW11X2xpc3QsIHNpZ21hPXNpZ21hX2xpc3QpCgojIENhbGN1bGF0ZSB0aGUgbGlrZWxpaG9vZCBvZiBldmVyeSBwb3NzaWJsZQojIGNvbWJpbmF0aW9uIGNhbmRpZGF0ZSBwYXJhbWV0ZXJzCnBvc3QkTEwgPC0gc2FwcGx5KDE6bnJvdyhwb3N0KSwgZnVuY3Rpb24oaSkgewogIHN1bShkbm9ybSgKICAgIGQyJGhlaWdodCwKICAgIG1lYW49cG9zdCRtdVtpXSwKICAgIHNkPXBvc3Qkc2lnbWFbaV0sCiAgICBsb2c9VFJVRSkpCiAgfSkKCiMgTXVsdGlwbHkgcHJpb3JzCnBvc3QkcHJvZCA8LSBwb3N0JExMICsgZG5vcm0ocG9zdCRtdSwgMTc4LCAyMCwgVFJVRSkgKwogIGR1bmlmKHBvc3Qkc2lnbWEsIDAsIDUwLCBUUlVFKQoKIyBDYWxjdWxhdGVkIHNjYWxlZCBwb3N0ZXJpb3IKcG9zdCRwcm9iIDwtIGV4cChwb3N0JHByb2QgLSBtYXgocG9zdCRwcm9kKSkKCiMgQ29udG91ciBwbG90IG9mIHByb2JhYmlsaXRpZXMKZ2dwbG90KHBvc3QsIGFlcyh4PW11LCB5PXNpZ21hLCB6PXByb2IpKSArIAogIGdlb21fY29udG91cihhZXMoY29sb3VyPXN0YXQobGV2ZWwpKSkgKyAKICBjb29yZF9jYXJ0ZXNpYW4oeGxpbT1jKDE0MCwgMTYwKSwgeWxpbT1jKDQsIDkpKQpgYGAKCk5vdyBsZXRzIHNhbXBsZSBmcm9tIHRoZSBwb3N0ZXJpb3IuCmBgYHtyfQojIFNhbXBsZSBjYW5kaWRhdGUgdmFsdWVzIGJ5IHBvc3RlcmlvcgpzYW1wbGVfcm93cyA8LSBzYW1wbGUoMTpucm93KHBvc3QpLCBzaXplPTFlNCwgcmVwbGFjZT1UUlVFLCBwcm9iPXBvc3QkcHJvYikKCiMgUGxvdCB0aGVpciAyZCBoaXN0b2dyYW0KZGF0YV9mcmFtZSgKICBzYW1wbGVfbXUgPSBwb3N0JG11W3NhbXBsZV9yb3dzXSwKICBzYW1wbGVfc2lnbWEgPSBwb3N0JHNpZ21hW3NhbXBsZV9yb3dzXQopICU+JSAKICBnZ3Bsb3QoYWVzKHNhbXBsZV9tdSwgc2FtcGxlX3NpZ21hKSkgKwogICAgZ2VvbV9iaW4yZChiaW5zPTMwKSArIAogICAgY29vcmRfY2FydGVzaWFuKHhsaW09YygxNDAsIDE2MCksIHlsaW09Yyg0LCA5KSkKYGBgCgpMZXRzIGFsc28gZXhhbWluZSB0aGUgbWFyZ2lucy4KYGBge3J9CnJlcXVpcmUoZ3JpZEV4dHJhKQoKcDEgPC0gcG9zdCRtdVtzYW1wbGVfcm93c10gJT4lIAogIHFwbG90ICsgZ2d0aXRsZSgibXUiKQpwMiA8LSBwb3N0JHNpZ21hW3NhbXBsZV9yb3dzXSAlPiUgCiAgcXBsb3QgKyBnZ3RpdGxlKCJzaWdtYSIpCgpncmlkLmFycmFuZ2UocDEsIHAyLCBuY29sPTIpCmBgYAoKYG1hcGAgZnJvbSB0aGUgcmV0aGlua2luZyBwYWNrYWdlIGFsbG93cyB1cyB0byBzb2x2ZSBmb3IgdGhlIHBvc3RlcmlvciB1c2luZyBxdWFkcmF0aWMgYXBwcm94aW1hdGlvbiBpbnN0ZWFkIG9mIGdyaWQgYXBwcm94aW1hdGlvbnMuIEl0IHNob3dzIHVzIGdhdXNzaWFuIGFwcHJveGltYXRpb25zIGZvciBlYWNoIHBhcmFtZXRlcidzIG1hcmdpbmFsIGRpc3RyaWJ1dGlvbi4KYGBge3J9CmZsaXN0IDwtIGFsaXN0KAogIGhlaWdodCB+IGRub3JtKG11LCBzaWdtYSksCiAgbXUgfiBkbm9ybSgxNzgsIDIwKSwKICBzaWdtYSB+IGR1bmlmKDAsIDUwKQopCm00LjEgPC0gcmV0aGlua2luZzo6bWFwKGZsaXN0LCBkYXRhPWQyKQpwcmVjaXMobTQuMSkKYGBgCgpWYXJpYW5jZSBjb3ZhcmlhbmNlIG1hdHJpeCB0ZWxscyB1cyBob3cgZWFjaCBwYXJhbWV0ZXIgcmVsYXRlcyB0byBldmVyeSBvdGhlciBwYXJhbWV0ZXIgaW4gdGhlIHBvc3RlcmlvciBkaXN0cmlidXRpb24uCmBgYHtyfQp2Y292KG00LjEpCiMgVGhpcyBjYW4gYmUgZmFjdG9yZWQgaW4gdG8gYSB2ZWN0b3Igb2YgdmFyaWFuY2VzCmRpYWcodmNvdihtNC4xKSkKIyBPciBhIGNvcnJlbGF0aW9uIG1hdHJpeApjb3YyY29yKHZjb3YobTQuMSkpCmBgYAoKTGV0cyBkcmF3IGZyb20gdGhlIHBvc3RlcmlvciBieSBzYW1wbGluZyB2ZWN0b3JzIG9mIHZhbHVlcyBmcm9tIGEgbXVsdGktZGltZW5zaW9uYWwgZ2F1c3NpYW4gZGlzdHJpYnV0aW9uLgpgYGB7cn0KcG9zdCA8LSBleHRyYWN0LnNhbXBsZXMobTQuMSwgbj0xZTQpCnByZWNpcyhwb3N0KQpgYGAKCiMjIyBBZGRpbmcgYSBwcmVkaWN0b3IgdmFyaWFibGUsIHdlaWdodAoKSXQgbG9va3MgbGlrZSB3ZWlnaHQgd2lsbCBiZSBhIGdvb2QgcHJlZGljdG9yIGZvciBoZWlnaHQgaW4gYWR1bHRzLgpgYGB7cn0KZ2dwbG90KGQyLCBhZXMoeD13ZWlnaHQsIHk9aGVpZ2h0KSkgKyAKICBnZW9tX3BvaW50KHNoYXBlPTEpICsKICBnZW9tX3Ntb290aChtZXRob2Q9bG0pCmBgYAoKVGhlIGdvYWwgaGVyZSBpcyB0byBtYWtlIHRoZSBwYXJhbWV0ZXIgZm9yIHRoZSBtZWFuIG9mIGEgR2F1c3NpYW4gZGlzdHJpYnV0aW9uLCAkXG11JCwgaW50byBhIGxpbmVhciBmdW5jdGlvbiBvZiB0aGUgcHJlZGljdG9yIHZhcmlhYmxlLiBUaGlzIGVuY29kZXMgdGhlIGFzc3VtcHRpb24gdGhhdCB0aGUgcHJlZGljdG9yIHZhcmlhYmxlIGhhcyBhIHBlcmZlY3RseSBjb25zdGFudCBhbmQgYWRkaXRpdmUgcmVsYXRpb25zaGlwIHRvIHRoZSBtZWFuIG9mIHRoZSBvdXRjb21lLgpgYGB7cn0KbTQuMyA8LSByZXRoaW5raW5nOjptYXAoCiAgYWxpc3QoCiAgICBoZWlnaHQgfiBkbm9ybShtdSwgc2lnbWEpLAogICAgbXUgPC0gYSArIGIqd2VpZ2h0LAogICAgYSB+IGRub3JtKDE1NiwgMTAwKSwKICAgIGIgfiBkbm9ybSgwLCAxMCksCiAgICBzaWdtYSB+IGR1bmlmKDAsIDUwKSksCiAgZGF0YSA9IGQyKQojIERpc3BsYXkgdGFibGUgb2YgZXN0aW1hdGVzCnByZWNpcyhtNC4zLCBjb3JyPVRSVUUpCmBgYAoKVGhlIHRhYmxlIGFib3ZlIHRlbGxzIHVzIHRoYXQgYSBwZXJzb24gMSBrZyBoZWF2aWVyIGlzIGV4cGVjdGVkIHRvIGJlIC45IGNtIHRhbGxlci4gSXQgYWxzbyB0ZWxscyB1cyB0aGF0IGEgcGVyc29uIHdpdGggMCB3ZWlnaHQgaXMgMTE0IGNtIHRhbGwsIHdoaWNoIGlzIG9mIGNvdXJzZSBub25zZW5zZTsgVGhpcyBpcyB3aHkgaXQgaXMgdXN1YWxseSBpbXBvcnRhbnQgdG8gaGF2ZSB2ZXJ5IHdlYWsgcHJpb3JzIGZvciBpbnRlcmNlcHRzIGluIG1hbnkgY2FzZXMuIEZpbmFsbHksIGl0IGV4cGxhaW5zIHRoYXQgOTUlIG9mIHBsYXVzaWJsZSBoZWlnaHRzIGxpZSB3aXRoaW4gMTAgY20gb2YgdGhlIG1lYW4uCgpDZW50ZXJpbmcgaXMgdGhlIHByb2NlZHVyZSBvZiBzdWJ0cmFjdGluZyB0aGUgbWVhbiBvZiBhIHZhcmlhYmxlIGZyb20gZWFjaCB2YWx1ZS4gCmBgYHtyfQpkMiR3ZWlnaHQuYyA8LSBkMiR3ZWlnaHQgLSBtZWFuKGQyJHdlaWdodCkKCm00LjQgPC0gcmV0aGlua2luZzo6bWFwKAogIGFsaXN0KAogICAgaGVpZ2h0IH4gZG5vcm0obXUsIHNpZ21hKSwKICAgIG11IDwtIGEgKyBiKndlaWdodC5jLAogICAgYSB+IGRub3JtKDE3OCwgMTAwKSwKICAgIGIgfiBkbm9ybSgwLCAxMCksCiAgICBzaWdtYSB+IGR1bmlmKDAsIDUwKSksCiAgZGF0YSA9IGQyKQojIERpc3BsYXkgdGFibGUgb2YgZXN0aW1hdGVzCnByZWNpcyhtNC40LCBjb3JyPVRSVUUpCmBgYAoKQnkgKipjZW50ZXJpbmcqKiBvdXIgcHJlZGljdG9yIHZhcmlhYmxlLCB0aGUgaW50ZXJwcmV0YXRpb24gb2YgdGhlIGludGVyY2VwdCBoYXMgbm93IGJlY29tZSB0aGUgZXhwZWN0ZWQgdmFsdWUgb2YgdGhlIG91dGNvbWUgd2hlbiB0aGUgcHJlZGljdG9yIGlzIGF0IGl0cyBhdmVyYWdlIHZhbHVlLiBUaGlzIG1ha2VzIGludGVycHJldGluZyB0aGUgaW50ZXJjZXB0IGEgbG90IGVhc2llci4gCgojIyMgVmlzdWFsaXppbmcgb3VyIGVzdGltYXRlcwoKV2UgY2FuIHBsb3Qgb3VyIG1heGltdW0gYSBwb3N0ZXJpb3JpIGZpdCBoZXJlCgpgYGB7cn0KcCA8LSBnZ3Bsb3QoZDIsIGFlcyh4PXdlaWdodCwgeT1oZWlnaHQpKSArCiAgZ2VvbV9wb2ludChzaGFwZT0xLCBzaXplPTIpICsKICBnZW9tX2FibGluZShhZXMoaW50ZXJjZXB0ID0gY29lZihtNC4zKVsiYSJdLAogICAgICAgICAgICAgICAgICBzbG9wZSA9IGNvZWYobTQuMylbImIiXSkpCnAKYGBgCgpUaGlzIGxpbmUgb25seSByZXByZXNlbnRzIG9uZSBwaWVjZSBvZiB0aGUgcG9zdGVyaW9yLiBXZSBjb3VsZCBpbnN0ZWFkIHNhbXBsZSB2YWx1ZXMgb2YgJFxhbHBoYSQgYW5kICRcYmV0YSQgZnJvbSB0aGUgcG9zdGVyaW9yIHRvIGdlbmVyYXRlIG1hbnkgbGluZXMgdG8gYnVpbGQgYm91bmRzIG9mIHVuY2VydGFpbnR5LiBXZSdsbCBkbyB0aGlzIHdpdGggYSBtdWNoIHNtYWxsZXIgZGF0YXNldCB0byBtYWduaWZ5IGVmZmVjdHMuCmBgYHtyfQpOIDwtIDEwCmROIDwtIGQyWyAxOk4gLCBdCm1OIDwtIHJldGhpbmtpbmc6Om1hcCgKICBhbGlzdCgKICAgIGhlaWdodCB+IGRub3JtKCBtdSAsIHNpZ21hICkgLAogICAgbXUgPC0gYSArIGIqd2VpZ2h0ICwKICAgIGEgfiBkbm9ybSggMTc4ICwgMTAwICkgLAogICAgYiB+IGRub3JtKCAwICwgMTAgKSAsCiAgICBzaWdtYSB+IGR1bmlmKCAwICwgNTAgKQopICwgZGF0YT1kTiApCgpwb3N0IDwtIGV4dHJhY3Quc2FtcGxlcyhtTiwgbj0yMCkKaGVhZChwb3N0KQpgYGAKCmBgYHtyfQpwIDwtIGdncGxvdChkTiwgYWVzKHdlaWdodCwgaGVpZ2h0KSkgKwogIGdlb21fcG9pbnQoc2hhcGU9MSkKCmZvciAoaSBpbiAxOjIwKSB7CiAgcCA8LSBwICsgZ2VvbV9hYmxpbmUoaW50ZXJjZXB0ID0gcG9zdCRhW2ldLCBzbG9wZSA9IHBvc3QkYltpXSwgYWxwaGEgPSBJKDMvMTApKQp9CgpwCmBgYAoKVG8gZG8gdGhpcyBtb3JlIGdlbmVyYWxseToKCjEuIFNhbXBsZSBwYXJhbWV0ZXIgdmFsdWVzIGZyb20gdGhlIHBvc3Rlcmlvci4KMi4gVXNlIHRoZSBwYXJhbWV0ZXIgdmFsdWVzIHRvIGNyZWF0ZSBwcmVkaWN0aW9ucyBvZiB0aGUgcGFyYW1ldGVycyBvZiBpbnRlcmVzdCwgaW5jbHVkaW5nIHRoZSBmaW5hbCBvdXRjb21lLgozLiBVc2Ugc3VtbWFyeSBmdW5jdGlvbnMgbGlrZSBtZWFuLCBIRE9QLCBQSSB0byBmaW5kIGF2ZXJhZ2VzIGFuZCBsb3dlciBhbmQgdXBwZXIgYm91bmRzIG9mIG91ciB1bmNlcnRhaW50eS4KCgpgYGB7cn0KIyBEcmF3IHBhcmFtZXRlciB2YWx1ZXMgZnJvbSB0aGUgcG9zdGVyaW9yIGRpc3RyaWJ1dGlvbgpwb3N0ZXJpb3Jfc2FtcGxlcyA8LSBleHRyYWN0LnNhbXBsZXMobTQuMykKCiMgQ3JlYXRlIG91dCB4LWF4aXMgb2YgaW50ZXJlc3QKd2VpZ2h0c19heGlzIDwtIDI1OjcwCgojIEZ1bmN0aW9uIHRvIGNhbGN1bGF0ZSBmaXQgdmFsdWVzIGZvciBcbXUKbXVfbGluayA8LSBmdW5jdGlvbih3ZWlnaHQpIHBvc3QkYSArIHBvc3QkYip3ZWlnaHQKbXUgPC0gcHVycnI6Om1hcCh3ZWlnaHRzX2F4aXMsIG11X2xpbmspCgojIEZ1bmN0aW9uIHRvIGNhbGN1bGF0ZSBmaXQgaGVpZ2h0cwojIFRoaXMgaW5jb3Jwb3JhdGVzIHZhcmlhYmlsaXR5IGFzIHdlbGwKc2ltdWxhdGUgPC0gZnVuY3Rpb24od2VpZ2h0LCBwYXJhbWV0ZXJfc2FtcGxlcykgewogIHJub3JtKG4gPSBucm93KHBhcmFtZXRlcl9zYW1wbGVzKSwKICAgICAgICBtZWFuID0gcGFyYW1ldGVyX3NhbXBsZXMkYSArIHBhcmFtZXRlcl9zYW1wbGVzJGIgKiB3ZWlnaHQsCiAgICAgICAgc2QgPSBwYXJhbWV0ZXJfc2FtcGxlcyRzaWdtYSl9CnNpbXVsYXRlZF9oZWlnaHRzIDwtIHB1cnJyOjptYXAod2VpZ2h0c19heGlzLCB+IHNpbXVsYXRlKC54LCBwb3N0ZXJpb3Jfc2FtcGxlcykpCgojIENvbWJpbmUgYWxsIHRoaXMgaW5mb3JtYXRpb24gdG8gcGxvdApkYXRhX2ZyYW1lKAogICMgWCBheGlzCiAgd2VpZ2h0ID0gd2VpZ2h0c19heGlzLAogICMgRXN0aW1hdGVzIGZvciBtZWFuIGhlaWdodAogIG11X21lYW4gPSBtYXBfZGJsKG11LCBtZWFuKSwKICBtdV9sb3dlciA9IG1hcF9kYmwobXUsIH5IUERJKC54LCAuODkpWzFdKSwKICBtdV91cHBlciA9IG1hcF9kYmwobXUsIH5IUERJKC54LCAuODkpWzJdKSwKICAjIEVzdGltYXRlcyBmb3IgaGVpZ2h0CiAgaGVpZ2h0X21lYW4gPSBtYXBfZGJsKHNpbXVsYXRlZF9oZWlnaHRzLCBtZWFuKSwKICBoZWlnaHRfbG93ID0gbWFwX2RibChzaW11bGF0ZWRfaGVpZ2h0cywgfkhQREkoLngsIC44OSlbMV0pLAogIGhlaWdodF9oaWdoID0gbWFwX2RibChzaW11bGF0ZWRfaGVpZ2h0cywgfkhQREkoLngsIC44OSlbMl0pKSAlPiUgCiAgZ2dwbG90KGRhdGE9LikgKwogIGdlb21fcG9pbnQoaW5oZXJpdC5hZXMgPSBGQUxTRSwgZGF0YSA9IGQyLAogICAgICAgICAgICAgYWVzKHg9d2VpZ2h0LCB5PWhlaWdodCksIHNoYXBlPTEpICsKICBnZW9tX2xpbmUoYWVzKHg9d2VpZ2h0LCB5PW11X21lYW4pKSArCiAgZ2VvbV9yaWJib24oYWVzKHggPSB3ZWlnaHQsIHltaW49bXVfbG93ZXIsIHltYXg9bXVfdXBwZXIpLCBhbHBoYT0uNSkgKwogIGdlb21fcmliYm9uKGFlcyh4ID0gd2VpZ2h0LCB5bWluID0gaGVpZ2h0X2xvdywgeW1heCA9IGhlaWdodF9oaWdoKSwgYWxwaGE9LjMpICsKICBnZ3RpdGxlKCJUaGUgZWZmZWN0IG9mIHdlaWdodCBvbiBoZWlnaHQiLCBzdWJ0aXRsZSA9ICJVc2luZyA4OSUgaGlnaGVzdCBwb3N0ZXJpb3IgZGVuc2l0eSBpbnRlcnZhbHMiKSArCiAgY29vcmRfY2FydGVzaWFuKHhsaW09YygzMCwgNjApLHlsaW09YygxMzUsIDE4MCkpCmBgYAoKIyMgUG9seW5vbWlhbCBSZWdyZXNzaW9uCgpTdGFuZGFyZGl6aW5nIG91ciBwcmVkaWN0b3JzIG1lYW5zIHNjYWxpbmcgYnV0IGFsc28gZGl2aWRpbmcgYnkgdGhlIHN0YW5kYXJkIGRldmlhdGlvbi4gVGhpcyBjaGFuZ2VzIG91ciBpbnRlcnByZXRhdGlvbiBvZiBjb250aW51b3VzIHZhcmlhYmxlcyB0byBtZWFuIHRoZSB1bml0IGNoYW5nZSBpbiBvdXRjb21lIGZyb20gYSBjYWhuZ2Ugb2Ygb25lIHN0YW5kYXJkIGRldmlhdGlvbi4gCgpgYGB7cn0KIyBTdGFuZGFyZGl6ZSBjb250aW51b3VzIHZhcmlhYmxlCnN0YW5kYXJkaXplIDwtIGZ1bmN0aW9uKHZlYykgewogIGNlbnRlcmVkIDwtIHZlYyAtIG1lYW4odmVjKQogIGNlbnRlcmVkIC8gc2QodmVjKQp9CmQkd2VpZ2h0LnMgPC0gc3RhbmRhcmRpemUoZCR3ZWlnaHQpCgojIERlZmluZSBhbm90aGVyIHBvbHlub21pYWwgdGVybQpkJHdlaWdodC5zMiA8LSBkJHdlaWdodC5zXjIKCm00LjUgPC0gbWFwKAogIGFsaXN0KAogICAgIyBPdXRjb21lCiAgICBoZWlnaHQgfiBkbm9ybShtdSwgc2lnbWEpLAogICAgIyBSZWdyZXNzaW9uIHNwZWNpZmljYXRpb24KICAgIG11IDwtIGEgKyBiMSp3ZWlnaHQucyArIGIyKndlaWdodC5zMiwKICAgIGEgfiBkbm9ybSgxNzggLCAxMDApLAogICAgIyBUaGVzZSBwcmlvcnMgYXJlIHZlcnkgd2VhayEKICAgIGIxIH4gZG5vcm0oMCAsIDEwKSwKICAgIGIyIH4gZG5vcm0oMCAsIDEwKSwKICAgIHNpZ21hIH4gZHVuaWYoMCAsIDUwKSksCiAgZGF0YT1kICkKcHJlY2lzKCBtNC41ICkKYGBgCgojIyBCLVNwbGluZXMKU3BsaW5lcyBhcmUgYW4gYWx0ZXJuYXRpdmUgdG8gbW9kZWxpbmcgbm9uLWxpbmVhcml0eS4gV2UnbGwgc2hvdyB0aGlzIHdpdGggc29tZSBjaGVycnkgYmxvc3NvbSBkYXRhLgpgYGB7cn0KbGlicmFyeShyZXRoaW5raW5nKQpsaWJyYXJ5KHJjYXJ0b2NvbG9yKQpjb2xvdXJfdGhlbWUgPC0gIkJ1cmdZbCIKcGFsZXR0ZSA8LSBjYXJ0b19wYWwoNywgY29sb3VyX3RoZW1lKQoKZGF0YShjaGVycnlfYmxvc3NvbXMpCmQgPC0gY2hlcnJ5X2Jsb3Nzb21zCgpkICU+JSAKICBzZWxlY3QoeWVhciwgdGVtcCkgJT4lIAogIGdncGxvdChhZXMoeWVhciwgdGVtcCkpICsKICAgIGdlb21fbGluZShjb2xvdXIgPSBwYWxldHRlWzddKSArCiAgICB0aGVtZV9idXJneWwoKQogICAgCmBgYAoKU3BsaW5lcyB3b3JrIGJ5IGRlZmluaW5nIGtub3QgcG9pbnRzIGFuZCBmaXR0aW5nIGxvY2FsIG1vZGVscyBiZXR3ZWVuIHRoZW0uIEEgY29tbW9uIHdheSB0byBwaWNrIHRoZXNlIGtub3QgcG9pbnRzIGlzIHRvIGV2ZW5seSBzcGFjZSB0aGVtIGFjcm9zcyB0aGUgcHJvYmFiaWxpdHkgZGVuc2l0eSBvZiB0aGUgdmFyaWFibGUuIENyb3NzIHZhbGlkYXRpb24gY2FuIHVsdGltYXRlbHkgYmUgdXNlZCB0byB2YWxpZGF0ZSB0aGUgbnVtYmVyIGFuZCBzcGFjaW5nLgpgYGB7cn0KZDIgPC0gZCAlPiUgCiAgZmlsdGVyKCFpcy5uYSh0ZW1wKSkKbnVtX2tub3RzIDwtIDE1Cmtub3RfbGlzdCA8LSBxdWFudGlsZShkMiR5ZWFyLCBwcm9icyA9IHNlcSgwLCAxLCBsZW5ndGgub3V0ID0gbnVtX2tub3RzKSkKZDIgJT4lIAogIGdncGxvdChhZXMoeWVhcikpICsKICBnZW9tX2RlbnNpdHkoY29sb3VyID0gInRyYW5zcGFyZW50IiwgZmlsbCA9IHBhbGV0dGVbNV0pICsKICBnZW9tX3BvaW50KGRhdGEgPSBkYXRhX2ZyYW1lKHk9IHJlcCgwLCBsZW5ndGgoa25vdF9saXN0KSksIHg9a25vdF9saXN0KSwKICAgICAgICAgICAgIG1hcHBpbmcgPSBhZXMoeCwgeSksIAogICAgICAgICAgICAgc2l6ZSA9IDMsIGNvbG91ciA9IHBhbGV0dGVbN10pICsKICB0aGVtZV9idXJneWwoKSArIGdndGl0bGUoIktub3QgcG9pbnQgYWxsb2NhdGlvbiIpCmBgYAoKTm93IHdlIHBpY2sgdGhlIHBvbHlub21pYWwgZGVncmVlLiBUaGlzIGRldGVybWluZXMgaG93IG1hbnkgYmFzaXMgcGFyYW1ldGVycyBhdCBvbmNlIGFyZSBpbnRlcmFjdGluZyB3aXRoIGEgZ2l2ZW4gcG9pbnQuIFRoZSBgc3BsaW5lc2AgbGlicmFyeSBpbiBSIGNhbiBoZWxwIHVzIHdpdGggdGhpcy4gVGhlIGNvZGUgYmVsb3cgd2lsbCBjcmVhdGUgMTcgY29sdW1ucyBpbiBwbGFjZSBvZiB5ZWFyIHRoYXQgcmVzdWx0IGZyb20gaXRzIGNoYW5nZSBvZiBiYXNpcy4KYGBge3J9CmxpYnJhcnkoc3BsaW5lcykKbGlicmFyeShicm1zKQoKIyBTcGxpdCB5ZWFyIGluIHRvIGJhc2lzIGZ1bmN0aW9uCkIgPC0gYnMoZDIkeWVhciwKICAgICAgICBrbm90cyA9IGtub3RfbGlzdFstYygxLCBudW1fa25vdHMpXSwKICAgICAgICBkZWdyZWUgPSAzLCBpbnRlcmNlcHQgPSBUUlVFKSAlPiUgCiAgYXNfZGF0YV9mcmFtZSAlPiUgCiAgc2V0X25hbWVzKC4sIHBhc3RlMCgieWVhcl8iLG5hbWVzKC4pKSkKCiMgQmluZCB0byBvcmlnaW5hbCBkYXRhIGZyYW1lCmQzIDwtIGQyICU+JSAKICBzZWxlY3QodGVtcCkgJT4lIAogIGJpbmRfY29scyhCKQpkMyAlPiUgZ2xpbXBzZQoKZml0LnNwbGluZSA8LSBicm0odGVtcCB+IC4sIGZhbWlseSA9IGdhdXNzaWFuKCksIGRhdGEgPSBkMywKICAgICAgICAgICAgICAgICAgcHJpb3IgPSBjKHByaW9yKG5vcm1hbCg2LCAxMCksIGNsYXNzID0gSW50ZXJjZXB0KSwKICAgICAgICAgICAgICAgICAgICAgICAgICAgIHByaW9yKG5vcm1hbCgwLCAxKSwgY2xhc3MgPSBiKSwKICAgICAgICAgICAgICAgICAgICAgICAgICAgIHByaW9yKGV4cG9uZW50aWFsKDEpLCBjbGFzcyA9IHNpZ21hKSksCiAgICAgICAgICAgICAgICAgIHJlZnJlc2ggPSAwKQpwb3N0ZXJpb3Jfc3VtbWFyeShmaXQuc3BsaW5lKQpgYGAKCkV4dHJhY3QgdGhlIHByZWRpY3RlZCBtZWFucyBhbmQgdGhlIGVycm9yIHRlcm0uCmBgYHtyfQptdSA9IGFzX2RhdGFfZnJhbWUoZml0dGVkKGZpdC5zcGxpbmUpKQp0ZW1wX2hhdCA9IGFzX2RhdGFfZnJhbWUocHJlZGljdChmaXQuc3BsaW5lKSkKCihkaW1lbnNpb25zIDwtIGxpc3QobXUgPSBkaW0obXUpLCAKICAgICAgICAgICAgICAgICAgIHRlbXBfaGF0ID0gZGltKHRlbXBfaGF0KSkpCmBgYAoKUGxvdCB0aGUgcG9zdGVyaW9yIHByZWRpY3Rpb25zLgpgYGB7cn0KZGF0YV9mcmFtZSgKICB0ZW1wID0gZmlsdGVyKGQsICFpcy5uYSh0ZW1wKSkkdGVtcCwKICB5ZWFyID0gZmlsdGVyKGQsICFpcy5uYSh0ZW1wKSkkeWVhciwKICB0ZW1wX211ID0gbXUkRXN0aW1hdGUsCiAgdGVtcF9tdV9sYiA9IG11JFEyLjUsCiAgdGVtcF9tdV91YiA9IG11JFE5Ny41LAogIHRlbXBfaGF0X3AgPSB0ZW1wX2hhdCRFc3RpbWF0ZSwKICB0ZW1wX2hhdF9sYiA9IHRlbXBfaGF0JFEyLjUsCiAgdGVtcF9oYXRfdWIgPSB0ZW1wX2hhdCRROTcuNSkgJT4lIAogIGdncGxvdChhZXMoeCA9IHllYXIsIHkgPSB0ZW1wKSkgKwogICAgZ2VvbV9qaXR0ZXIoYWxwaGEgPSAuNCwgY29sb3VyID0gcGFsZXR0ZVs2XSkgKwogICAgZ2VvbV9yaWJib24oYWVzKHltaW4gPSB0ZW1wX211X2xiLCB5bWF4ID0gdGVtcF9tdV91YiksIGZpbGwgPSBhbHBoYShwYWxldHRlWzRdLCAuNykpICsKICAgIGdlb21fcmliYm9uKGFlcyh5bWluID0gdGVtcF9oYXRfbGIsIHltYXggPSB0ZW1wX2hhdF91YiksIGZpbGwgPSBhbHBoYShwYWxldHRlWzRdLCAuNCkpICsKICAgIGdlb21fbGluZShhZXMoeSA9IHRlbXBfaGF0X3ApLCBjb2xvdXIgPSBwYWxldHRlWzddKSArCiAgICB0aGVtZV9idXJneWwoKSArCiAgICBnZ3RpdGxlKCJQb3N0ZXJpb3IgcHJlZGljdGlvbiBjaGVjayIsCiAgICAgICAgICAgICJGaXR0aW5nIHRlbXBlcmF0dXJlIHRvIGEgMTUga25vdHRlZCBjdWJpYyBzcGxpbmUgb2YgeWVhciIpCmBgYAoKQWN0dWFsbHksIGJybXMgY2FuIHVzaW5nIHRoZSBgc2Agb3IgYHQyYCBmdW5jdGlvbnMgZnJvbSB0aGUgYG1nY3ZgIHBhY2thZ2UgdG8gYWNjb21wbGlzaCB0aGUgc2FtZSB0aGluZy4KYGBge3J9CmZpdC5zcGxpbmUuMiA8LSBicm0odGVtcCB+IHMoeWVhciwgayA9IDE1KSwgIyAxNSBrbm90cwogICAgICAgICAgICAgICAgICAgIGRhdGEgPSBkMiwgZmFtaWx5ID0gZ2F1c3NpYW4oKSwKICAgICAgICAgICAgICAgICAgICBwcmlvciA9IGMocHJpb3Iobm9ybWFsKDYsIDEwKSwgY2xhc3MgPSBJbnRlcmNlcHQpLAogICAgICAgICAgICAgICAgICAgICAgICAgICAgcHJpb3Iobm9ybWFsKDAsIDEpLCBjbGFzcyA9IGIpLAogICAgICAgICAgICAgICAgICAgICAgICAgICAgcHJpb3IoZXhwb25lbnRpYWwoMSksIGNsYXNzID0gc2lnbWEpKSkKCm11ID0gYXNfZGF0YV9mcmFtZShmaXR0ZWQoZml0LnNwbGluZS4yKSkKdGVtcF9oYXQgPSBhc19kYXRhX2ZyYW1lKHByZWRpY3QoZml0LnNwbGluZS4yKSkKCmRhdGFfZnJhbWUoCiAgdGVtcCA9IGZpbHRlcihkLCAhaXMubmEodGVtcCkpJHRlbXAsCiAgeWVhciA9IGZpbHRlcihkLCAhaXMubmEodGVtcCkpJHllYXIsCiAgdGVtcF9tdSA9IG11JEVzdGltYXRlLAogIHRlbXBfbXVfbGIgPSBtdSRRMi41LAogIHRlbXBfbXVfdWIgPSBtdSRROTcuNSwKICB0ZW1wX2hhdF9wID0gdGVtcF9oYXQkRXN0aW1hdGUsCiAgdGVtcF9oYXRfbGIgPSB0ZW1wX2hhdCRRMi41LAogIHRlbXBfaGF0X3ViID0gdGVtcF9oYXQkUTk3LjUpICU+JSAKICBnZ3Bsb3QoYWVzKHggPSB5ZWFyLCB5ID0gdGVtcCkpICsKICAgIGdlb21faml0dGVyKGFscGhhID0gLjQsIGNvbG91ciA9IHBhbGV0dGVbNl0pICsKICAgIGdlb21fcmliYm9uKGFlcyh5bWluID0gdGVtcF9tdV9sYiwgeW1heCA9IHRlbXBfbXVfdWIpLCBmaWxsID0gYWxwaGEocGFsZXR0ZVs0XSwgLjcpKSArCiAgICBnZW9tX3JpYmJvbihhZXMoeW1pbiA9IHRlbXBfaGF0X2xiLCB5bWF4ID0gdGVtcF9oYXRfdWIpLCBmaWxsID0gYWxwaGEocGFsZXR0ZVs0XSwgLjQpKSArCiAgICBnZW9tX2xpbmUoYWVzKHkgPSB0ZW1wX2hhdF9wKSwgY29sb3VyID0gcGFsZXR0ZVs3XSkgKwogICAgdGhlbWVfYnVyZ3lsKCkgKwogICAgZ2d0aXRsZSgiUG9zdGVyaW9yIHByZWRpY3Rpb24gY2hlY2siLAogICAgICAgICAgICAiRml0dGluZyB0ZW1wZXJhdHVyZSB0byBhIDE1IGtub3R0ZWQgY3ViaWMgc3BsaW5lIG9mIHllYXIiKQoKYGBgCgoKIyMgT24gdHJhbnNmb3JtYXRpb25zCgpDZW50ZXJpbmcgYW5kIHNjYWxpbmcgYSBjb250aW51b3VzIHZhcmlhYmxlIGluY3JlYXNlcyBpbnRlcnByZXRhYmlsaXR5LiBEZW1lYW5pbmcgYSBwcmVkaWN0b3IgdmFyaWFibGUgc3VidHJhY3RzIGV2ZXJ5IHZhbHVlIGJ5IGl0J3MgbWVhbi4gQnkgZGVtZWFuaW5nIGEgcHJlZGljdG9yIHZhcmlhYmxlcywgdGhlIGludGVyY2VwdCByZXByZXNlbnRzIHRoZSB2YWx1ZSBvZiB0aGUgb3V0Y29tZSB3aGVuIGFsbCB2YXJpYWJsZXMgYXJlIGF0IHRoZWlyIG1lYW4uIE9uZSBzdGFuZGFyZGl6ZXMgYSB2YXJpYWJsZSBieSBkaXZpZGluZyBldmVyeSB2YWx1ZSBieSB0aGUgdmFyaWFibGVzIHN0YW5kYXJkIGRldmlhdGlvbi4gU3RhbmRhcmRpemluZyBjYXVzZXMgdGhlIGludGVycHJldGF0aW9uIG9mIGNvZWZmaWNpZW50cyB0byBiZSBpbiB0ZXJtcyBvZiB1bml0IGNoYW5nZXMgaW4gU0Qgb2YgdGhlIHZhcmlhYmxlLgoKCkxvZyB0cmFuc2Zvcm1pbmcgdGhlIG91dGNvbWUgdmFyaWFibGUgbWFrZXMgeSBvdXIgY29lZmZpY2llbnRzIHJlcHJlc2VudCBwZXJjZW50IGluY3JlYXNlcyBpbiB0aGUgb3V0Y29tZS4gTG9nIHRyYW5zZm9ybWluZyBhIHByZWRpY3RvciB2YXJpYWJsZSBkZXNjcmliZXMgaG93IGEgc2luZ2xlIHBlcmNlbnRhZ2UgcG9pbnQgaW5jcmVhc2UgaW4gdGhlIHByZWRpY3RvciBhZmZlY3RzIHkuIExvZyB0cmFuc2Zvcm1pbmcgYm90aCB0aGUgcHJlZGljdG9yIGFuZCBvdXRjb21lIGRlc2NyaWJlcyBob3cgcGVyY2VudGFnZSBjaGFuZ2VzIGluIG9uZSBjcmVhdGUgcGVyY2VudGFnZSBjaGFuZ2VzIGluIHRoZSBvdGhlci4KCgojIEhvbWV3b3JrCgojIyBNZWRpdW0KIyMjIDRNMQpGb3IgdGhlIG1vZGVsIGRlZmluaXRpb24gYmVsb3csIHNpbXVsYXRlIG9ic2VydmVkIGhlaWdodHMgZnJvbSB0aGUgcHJpb3IuCmBgYHtyfQpudW1fb2JzZXJ2YXRpb25zIDwtIDFlNAoKcHJpb3JfbWVhbiA8LSBybm9ybShudW1fb2JzZXJ2YXRpb25zLCAwLCAxMCkKcHJpb3Jfc2QgPC0gcnVuaWYobnVtX29ic2VydmF0aW9ucywgMCwgMTApCgpwcmlvcl9oZWlnaHRzIDwtIHJub3JtKG51bV9vYnNlcnZhdGlvbnMsIHByaW9yX21lYW4sIHByaW9yX3NkKQoKcXBsb3QocHJpb3JfaGVpZ2h0cykKYGBgCgoKIyMjIDRNMgpgYGB7cn0KbW9kZWwgPC0gYWxpc3QoCiAgeSB+IGRub3JtKG11LCBzaWdtYSksCiAgbXUgfiBkbm9ybSgwLCAxMCksCiAgc2lnbWEgfiBkdW5pZigwLCAxMCkKKQpgYGAKCgoKIyMgSGFyZAojIyMgNEgxCmBgYHtyfQpkYXRhKCJIb3dlbGwxIikKZCA8LSBIb3dlbGwxCmQkd2VpZ2h0IDwtIHN0YW5kYXJkaXplKGQkd2VpZ2h0KQpkCmBgYAoKCmBgYHtyfQptb2RlbCA8LSBtYXAoCiAgYWxpc3QoCiAgICBoZWlnaHQgfiBkbm9ybShtdSwgc2lnbWEpLAogICAgbXUgPC0gYSArIGIgKiB3ZWlnaHQsCiAgICBhIH4gZG5vcm0obWVhbiA9IDAsIHNkID0gMTAwKSwKICAgIGIgfiBkbm9ybShtZWFuID0gMCwgc2QgPSAxMCksCiAgICBzaWdtYSB+IGR1bmlmKG1pbiA9IDAsIG1heCA9IDY0KQogICksCiAgZGF0YSA9IGQKKQptb2RlbApgYGAKCmBgYHtyfQojIFN0YW5kYXJkaXplIG91ciBpbnB1dApuZXdfd2VpZ2h0cyA8LSBzdGFuZGFyZGl6ZShjKDQ2Ljk1LCA0My43MiwgNjQuNzgsIDMyLjU5LCA1NC42MykpCgojIERyYXcgcGFyYW1ldGVyIHZhbHVlcyBmcm9tIHRoZSBwb3N0ZXJpb3IgZGlzdHJpYnV0aW9uCnBvc3Rlcmlvcl9zYW1wbGVzIDwtIGV4dHJhY3Quc2FtcGxlcyhtb2RlbCkKCiMgRnVuY3Rpb24gdG8gY2FsY3VsYXRlIGZpdCBoZWlnaHRzCnNpbXVsYXRlX2hlaWdodHMgPC0gZnVuY3Rpb24od2VpZ2h0LCBwYXJhbWV0ZXJfc2FtcGxlcykgewogIHJub3JtKAogICAgbiA9IG5yb3cocGFyYW1ldGVyX3NhbXBsZXMpLAogICAgbWVhbiA9IHBhcmFtZXRlcl9zYW1wbGVzJGEgKyBwYXJhbWV0ZXJfc2FtcGxlcyRiICogd2VpZ2h0LAogICAgc2QgPSBwYXJhbWV0ZXJfc2FtcGxlcyRzaWdtYSl9CgojIEdlbmVyYXRlIGEgcHJlZGljdGlvbiBmb3IgZWFjaCBuZXcgd2VpZ2h0CiMgVXNlIGV2ZXJ5IHBhcmFtZXRlciBzZXQgaW4gdGhlIHBvc3RlcmlvciB0byBkbyBzbwpzaW11bGF0ZWRfaGVpZ2h0cyA8LSBwdXJycjo6bWFwKG5ld193ZWlnaHRzLCB+IHNpbXVsYXRlKC54LCBwb3N0ZXJpb3Jfc2FtcGxlcykpCgpkYXRhX2ZyYW1lKAogIGluZGl2aWR1YWwgPSAxOjUsCiAgd2VpZ2h0ID0gYyg0Ni45NSwgNDMuNzIsIDY0Ljc4LCAzMi41OSwgNTQuNjMpLAogIGV4cGVjdGVkX2hlaWdodCA9IG1hcF9kYmwoc2ltdWxhdGVkX2hlaWdodHMsIG1lYW4pLAogIGxvd2VyX3BpID0gbWFwX2RibChzaW11bGF0ZWRfaGVpZ2h0cywgflBJKC54KVsxXSksCiAgdXBwZXJfcGkgPSBtYXBfZGJsKHNpbXVsYXRlZF9oZWlnaHRzLCB+UEkoLngpWzJdKQopCmBgYAoKIyA0SDIKYGBge3J9CmRhdGEoIkhvd2VsbDEiKQpkIDwtIEhvd2VsbDEKZCAlPiUgCiAgZmlsdGVyKGFnZSA8IDE4KSAlPiUgCiAgZ2dwbG90KGFlcyh4ID0gd2VpZ2h0LCB5ID0gaGVpZ2h0KSkgKwogIGdlb21fcG9pbnQoKSArIGdlb21fc21vb3RoKG1ldGhvZD1sbSkgKwogIGdlb21fc21vb3RoKGNvbG91ciA9ICJvcmFuZ2UiKQpgYGAKCmBgYHtyfQpkICU+JSAKICBmaWx0ZXIoYWdlIDwgMTgpICU+JSAKICBtYXAoYWxpc3QoCiAgICBoZWlnaHQgfiBkbm9ybShtdSwgc2lnbWEpLAogICAgbXUgPC0gYSArIGIgKiB3ZWlnaHQsCiAgICBhIH4gZG5vcm0obWVhbiA9IDEwMCwgc2QgPSAxMDApLAogICAgYiB+IGRub3JtKG1lYW4gPSAwLCBzZCA9IDEwKSwKICAgIHNpZ21hIH4gZHVuaWYobWluID0gMCwgbWF4ID0gNTApCiAgKSwKICBkYXRhID0gLikgLT4KbW9kZWwKCnByZWNpcyhtb2RlbCkKYGBgCgpBIG9uZSB1bml0IGluY3JlYXNlIGluIHdlaWdodCBpbmNyZWFzZXMgaGVpZ2h0IGJ5IDIuNzIsIHRoZXJlZm9yZSBhIDEwIHVuaXQgaW5jcmVhc2Ugd2lsbCBtb3ZlIGhlaWdodCBieSAyNy4yLgoKYGBge3J9CiMgU2FtcGxlIHRoZSBwb3NldGVyaW9yCnBvc3Rlcmlvcl9zYW1wbGVzIDwtIGV4dHJhY3Quc2FtcGxlcyhtb2RlbCkKCiMgQ3JlYXRlIHgtYXhpcwpuZXdfd2VpZ2h0cyA8LSBzZXEoZnJvbSA9IDQsIHRvID0gNDUsIGxlbmd0aC5vdXQgPSAxMDApCgojIENhbGN1bGF0ZSB0aGUgbWVhbgpmbl9tZWFuIDwtIGZ1bmN0aW9uKHdlaWdodCkge3Bvc3Rlcmlvcl9zYW1wbGVzJGEgKyBwb3N0ZXJpb3Jfc2FtcGxlcyRiICogd2VpZ2h0fQpmaXRfbWVhbnMgPC0gcHVycnI6Om1hcChuZXdfd2VpZ2h0cywgZm5fbWVhbikKCiMgQ2FsY3VsYXRlIGhlaWdodHMKZm5faGVpZ2h0IDwtIGZ1bmN0aW9uKHdlaWdodCkgewogIHJub3JtKAogICAgbiA9IG5yb3cocG9zdGVyaW9yX3NhbXBsZXMpLAogICAgbWVhbiA9IHBvc3Rlcmlvcl9zYW1wbGVzJGEgKyBwb3N0ZXJpb3Jfc2FtcGxlcyRiICogd2VpZ2h0LAogICAgc2QgPSBwb3N0ZXJpb3Jfc2FtcGxlcyRzaWdtYSkKfQpmaXRfaGVpZ2h0cyA8LSBwdXJycjo6bWFwKG5ld193ZWlnaHRzLCBmbl9oZWlnaHQpCgpkYXRhX2ZyYW1lKAogIHdlaWdodHMgPSBuZXdfd2VpZ2h0cywKICBtZWFucyA9IG1hcF9kYmwoZml0X21lYW5zLCBtZWFuKSwKICBtZWFuc19sb3cgPSBtYXBfZGJsKGZpdF9tZWFucywgfiBIUERJKC54LCAuODkpWzFdKSwKICBtZWFuc19oaWdoID0gbWFwX2RibChmaXRfbWVhbnMsIH4gSFBESSgueCwgLjg5KVsyXSksCiAgaGVpZ2h0cyA9IG1hcF9kYmwoZml0X2hlaWdodHMsIG1lYW4pLAogIGhlaWdodHNfbG93ID0gbWFwX2RibChmaXRfaGVpZ2h0cywgfiBIUERJKC54LCAuODkpWzFdKSwKICBoZWlnaHRzX2hpZ2ggPSBtYXBfZGJsKGZpdF9oZWlnaHRzLCB+IEhQREkoLngsIC44OSlbMl0pLAopICU+JSAKZ2dwbG90KGFlcyh4ID0gd2VpZ2h0cykpICsKICBnZW9tX2xpbmUoYWVzKHkgPSBtZWFucykpICsKICBnZW9tX3JpYmJvbihhZXMoeW1pbiA9IG1lYW5zX2xvdywgeW1heCA9IG1lYW5zX2hpZ2gpLCBhbHBoYSA9IC40KSArCiAgZ2VvbV9yaWJib24oYWVzKHltaW4gPSBoZWlnaHRzX2xvdywgeW1heCA9IGhlaWdodHNfaGlnaCksIGFscGhhID0gLjIpICsKICBnZW9tX3BvaW50KGluaGVyaXQuYWVzID0gRkFMU0UsIGRhdGEgPSBmaWx0ZXIoZCwgYWdlIDwgMTgpLAogICAgICAgICAgICAgYWVzKHg9d2VpZ2h0LCB5ID0gaGVpZ2h0KSwgc2hhcGUgPSAxKSArCiAgeWxhYigiaGVpZ2h0IikgKyB4bGFiKCJ3ZWlnaHQiKSArIHRoZW1lX2J3KCkgKwogIHRoZW1lKHBhbmVsLmdyaWQubWFqb3IgPSBlbGVtZW50X2JsYW5rKCksIHBhbmVsLmdyaWQubWlub3IgPSBlbGVtZW50X2JsYW5rKCkpCmBgYAoKIyBCUk1TCgpgYGB7cn0KZGF0YShIb3dlbGwxKQpkIDwtIEhvd2VsbDEKZDIgPC0gZCAlPiUgCiAgZmlsdGVyKGFnZSA+PSAxOCkKZDIgJT4lIAogIHNraW0KYGBgCgoKYGBge3J9CmxpYnJhcnkoYnJtcykKCmI0LjEgPC0KICBicm0oZGF0YSA9IGQyLCBmYW1pbHkgPSBnYXVzc2lhbiwKICAgICAgaGVpZ2h0IH4gMSwKICAgICAgcHJpb3IgPSBjKHByaW9yKG5vcm1hbCgxNzgsIDIwKSwgY2xhc3MgPSBJbnRlcmNlcHQpLAogICAgICAgICAgICAgICAgcHJpb3IoY2F1Y2h5KDAsIDEpLCBjbGFzcyA9IHNpZ21hKSksCiAgICAgIGl0ZXIgPSAzMTAwMCwgd2FybXVwID0gMzAwMDAsIGNoYWlucyA9IDQsIGNvcmVzID0gNCkKcGxvdChiNC4xKQpgYGAKCkV4dHJhY3QgdmFyaWFuY2UgY292YXJpYW5jZSBtYXRyaXguCmBgYHtyfQojIEV4dHJhY3Qgc2FtcGxlcyBmcm9tIHRoZSBwb3N0ZXJpb3IKcG9zdCA8LSBwb3N0ZXJpb3Jfc2FtcGxlcyhiNC4xKQojIENvdmFyaWFuY2UgbWF0cml4CmNvdihwb3N0WywgMToyXSkKYGBgCgpMZXRzIGV4dHJhY3QgdGhlIEhEUEkKYGBge3J9CmxpYnJhcnkodGlkeWJheWVzKQoKcG9zdCAlPiUgCiAgc2VsZWN0KC1scF9fKSAlPiUgCiAgZ2F0aGVyKHBhcmFtZXRlcikgJT4lIAogIGdyb3VwX2J5KHBhcmFtZXRlcikgJT4lIAogIG1lZGlhbl9oZGkodmFsdWUpCmBgYAoKTm93IGxldHMgYWRkIGEgcHJlZGljdG9yLgpgYGB7ciBtZXNzYWdlPUZBTFNFLCB3YXJuaW5nPUZBTFNFfQpkMyA8LSBkMiAlPiUgCiAgbXV0YXRlKHdlaWdodCA9ICh3ZWlnaHQgLSBtZWFuKHdlaWdodCkpL3NkKHdlaWdodCkpICU+JSAKICBzZWxlY3QoaGVpZ2h0LCB3ZWlnaHQpCgpiNC4zIDwtIAogIGJybShkYXRhID0gZDMsIGZhbWlseSA9IGdhdXNzaWFuLAogICAgICBoZWlnaHQgfiAxICsgd2VpZ2h0LAogICAgICBwcmlvciA9IGMocHJpb3Iobm9ybWFsKDE1NiwgMTAwKSwgY2xhc3MgPSBJbnRlcmNlcHQpLAogICAgICAgICAgICAgICAgcHJpb3Iobm9ybWFsKDAsIDEwKSwgY2xhc3MgPSBiKSwKICAgICAgICAgICAgICAgIHByaW9yKHVuaWZvcm0oMCwgNTApLCBjbGFzcyA9IHNpZ21hKSksCiAgICAgIGl0ZXIgPSA0NjAwMCwgd2FybXVwID0gNDUwMDAsIGNoYWlucyA9IDQsIGNvcmVzID0gNCwKICAgICAgY29udHJvbCA9IGxpc3QoYWRhcHRfZGVsdGEgPSAwLjgsIAogICAgICAgICAgICAgICAgICAgICBtYXhfdHJlZWRlcHRoID0gMTApKQoKcGxvdChiNC4zKQpgYGAKCkV4YW1pbmUgY292YXJpYW5jZQpgYGB7cn0KcGFpcnMoYjQuMykKYGBgCgpDcmVhdGluZyBwcmVkaWN0aW9ucwpgYGB7cn0KcG9zdCA8LSBwb3N0ZXJpb3Jfc2FtcGxlcyhiNC4zKQoKbXVfYXRfNTAgPC0KICBwb3N0ICU+JSAKICB0cmFuc211dGUobXVfYXRfNTAgPSBiX0ludGVyY2VwdCArIGJfd2VpZ2h0ICogNTApCgptdV9hdF81MCAlPiUKICBnZ3Bsb3QoYWVzKHggPSBtdV9hdF81MCkpICsKICBnZW9tX2RlbnNpdHkoc2l6ZSA9IDAsIGZpbGwgPSAiZ3JleTc1IikgKwogIHN0YXRfcG9pbnRpbnRlcnZhbGgoYWVzKHkgPSAwKSwgCiAgICAgICAgICAgICAgICAgICAgICBwb2ludF9pbnRlcnZhbCA9IG1vZGVfaGRpLCAud2lkdGggPSAuOTUpICsKICBzY2FsZV95X2NvbnRpbnVvdXMoTlVMTCwgYnJlYWtzID0gTlVMTCkgKwogIGxhYnMoeCA9IGV4cHJlc3Npb24obXVbImhlaWdodCB8IHdlaWdodCA9IDUwIl0pKSArCiAgdGhlbWVfY2xhc3NpYygpCgpgYGAKCk5vdyB3ZSBleHRyYWN0IGNvbXBsZXRlIHByZWRpY3Rpb25zIGluY2x1ZGluZyB0aGUgbWVhbiBhbmQgYWRkZWQgdmFyaWFuY2UuCmBgYHtyfQp3ZWlnaHRfc2VxIDwtIHRpYmJsZSh3ZWlnaHQgPSBzZXEoZnJvbSA9IDI1LCB0byA9IDcwLCBieSA9IDEpKQoKcHJlZF9oZWlnaHQgPC0KICBwcmVkaWN0KGI0LjMsCiAgICAgICAgICBuZXdkYXRhID0gd2VpZ2h0X3NlcSkgJT4lCiAgYXNfdGliYmxlKCkgJT4lCiAgYmluZF9jb2xzKHdlaWdodF9zZXEpCiAgCnByZWRfaGVpZ2h0ICU+JQogIHNsaWNlKDE6NikKYGBgCgpgYGB7cn0KZ2dwbG90KHByZWRfaGVpZ2h0LCBhZXMod2VpZ2h0LCBFc3RpbWF0ZSkpICsKICBnZW9tX3JpYmJvbihhZXMoeW1pbiA9IFEyLjUsIHltYXggPSBROTcuNSksIGZpbGwgPSAiZ3JleTc1IikgKwogIGdlb21fbGluZSgpICsKICB0aGVtZShwYW5lbC5ncmlkID0gZWxlbWVudF9ibGFuaygpKQpgYGAKCldlIGNhbiBvdmVybGF5IHRoaXMgd2l0aCB0aGUgbWVhbiBpbnRlcnZhbCB0b28KYGBge3J9Cm11X3N1bW1hcnkgPC0KICBmaXR0ZWQoYjQuMywgCiAgICAgICAgIG5ld2RhdGEgPSB3ZWlnaHRfc2VxKSAlPiUKICBhc190aWJibGUoKSAlPiUKICBiaW5kX2NvbHMod2VpZ2h0X3NlcSkKCmdncGxvdChtdV9zdW1tYXJ5LCBhZXMod2VpZ2h0LCBFc3RpbWF0ZSkpICsKICBnZW9tX3JpYmJvbihhZXMoeW1pbiA9IFEyLjUsIHltYXggPSBROTcuNSksIGZpbGwgPSAiZ3JleTc1IikgKwogIGdlb21fbGluZSgpICsKICB0aGVtZShwYW5lbC5ncmlkID0gZWxlbWVudF9ibGFuaygpKQpgYGAKCmBgYHtyfQppbm5lcl9qb2luKHByZWRfaGVpZ2h0LCBtdV9zdW1tYXJ5LCBieSA9ICJ3ZWlnaHQiLCBzdWZmaXg9YygiX2hlaWdodCIsICJfbWVhbiIpKSAlPiUgc3RyCmBgYAoKYGBge3J9CmlubmVyX2pvaW4ocHJlZF9oZWlnaHQsIG11X3N1bW1hcnksIGJ5ID0gIndlaWdodCIsIHN1ZmZpeD1jKCJfaGVpZ2h0IiwgIl9tZWFuIikpICU+JSAKICBnZ3Bsb3QoYWVzKHggPSB3ZWlnaHQsIHkgPSBFc3RpbWF0ZV9oZWlnaHQpKSArCiAgZ2VvbV9yaWJib24oYWVzKHltaW4gPSBRMi41X2hlaWdodCwgeW1heCA9IFE5Ny41X2hlaWdodCksIGZpbGwgPSAiZ3JleTI1IikgKwogIGdlb21fcmliYm9uKGFlcyh5bWluID0gUTIuNV9tZWFuLCB5bWF4ID0gUTk3LjVfbWVhbiksIGZpbGwgPSAiZ3JleTQ1IikgKwogIGdlb21fbGluZSgpICsKICB0aGVtZShwYW5lbC5ncmlkID0gZWxlbWVudF9ibGFuaygpKQpgYGAKCgoKCgoKCgoKCgo=

Chapter 4 - Linear Models